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ABSTRACT 

We study of the collapse of a magnetized spherical star to a black hole in general relativity theory. 
The matter and gravitational fields are described by the exact Oppenheimer-Snyder solution for the 
collapse of a spherical, homogeneous dust ball. We adopt a “dynamical Cowling approximation” whereby 
the matter and the geometry (metric), while highly dynamical, are unaffected by the electromagnetic 
fields. The matter is assumed to be perfectly conducting and threaded by a dipole magnetic field 
at the onset of collapse. We determine the subsequent evolution of the magnetic and electric fields 
without approximation; the fields are determined analytically in the matter interior and numerically 
in the vacuum exterior. We apply junction conditions to match the electromagnetic fields across the 
stellar surface. We use this model to experiment with several coordinate gauge choices for handling 
spacetime evolution characterized by the formation of a black hole and the associated appearance of 
singularities. These choices range from “singularity avoiding” time coordinates to “horizon penetrating” 
time coordinates accompanied by black hole excision. The later choice enables us to integrate the 
electromagnetic fields arbitrarily far into the future. At late times the longitudinal magnetic field in 
the exterior has been transformed into a transverse electromagnetic wave; part of the electromagnetic 
radiation is captured by the hole and the rest propagates outward to large distances. The solution we 
present for our simple scenario can be used to test codes designed to treat more general evolutions of 
relativistic MHD fluids flowing in strong gravitational fields in dynamical spacetimes. 


1. INTRODUCTION 

Magnetic fields play a crucial role in determining the 
evolution of many relativistic objects. In a compan¬ 
ion paper (Baumgarte & Shapiro, 2002, hereafter Paper 
I) we assembled the complete set of MaxwelDEinstein- 
magnetohydrodynamic (MHD) equations that determine 
the self-consistent evolution of a relativistic, ideal MHD 
gas in a dynamical spacetime. Our goal was to set down a 
formulation of the equations that is suitable for numerical 
integration in full 3-1-1 dimensions. 

In this paper (hereafter Paper H) we use these gen¬ 
eral relativistic MHD equations to follow the gravitational 
collapse of a magnetized star to a black hole. We solve 
this problem in a “dynamical Cowling approximation”, 
whereby the matter and gravitational field are given by 
the Oppenheimer-Sndyer solution (Oppenheimer & Sny¬ 
der 1939) for the collapse from rest of a spherical dust ball 
to a Schwarzschild black hole. We assume that initially 
the star is threaded by a dipole magnetic field and deter¬ 
mine the subsequent evolution of the magnetic and electric 
fields with time, both in the matter interior and vacuum 
exterior. 

We are motivated to tackle this problem partly to accli¬ 
mate ourselves to the task of numerically integrating the 
coupled Maxwell-Einstein-relativistic MHD equations for 
a strong-held, dynamical spacetime. Though the back¬ 
ground geometry is entirely determined in this example, 
the spacetime is highly dynamical and characterized by 
a very strong (black hole) gravitational held. The goal 
is to solve for the electromagnetic helds to arbitrary late 
times everywhere in space. The restricted problem we pose 


here is sufficiently rich to allow a comparison of coordi¬ 
nate gauge choices for handling collapse to a black hole 
and dealing with the accompanying spacetime singularity. 
Such gauge choices range from “singularity avoiding” time 
coordinates to horizon penetrating time slices which allow 
for black hole excision. We experiment with these different 
choices below. 

The collapse of a magnetized star to a black hole is an 
astrophysically important problem that has been discussed 
previously. Some of the earliest treatments were based on 
a quasi-static approach (Ginzburg & Ozernoy 1965; An¬ 
derson & Gohen 1970; Zel’dovich & Novikov 1971) which 
was demonstrated to give the incorrect asymptotic decay 
of the fields with time (Price 1972a,b; see also Thorne 
1971). The general problem of magnetic star collapse to a 
black role obviously requires a detailed numerical treat¬ 
ment for solution of the coupled Maxwell - Einstein - 
relativistic MHD system. Only now is such a treatment 
feasible. In fact, we are unaware of any previous analysis 
that determines the complete evolution of both the interior 
and exterior electromagnetic fields during stellar collapse 
to a black hole, even for the restricted spherical collapse 
problem posed here. [In a pioneering paper, Wilson 1975 
used relativistic MHD to follow the collapse of the interior 
of a magnetized star, but treated the gravitational field 
as quasi-stationary]. Our solution should serve as a useful 
testbed for designing codes capable of handling more com¬ 
plicated scenarios involving relativistic MHD gravitational 
collapse. Such scenarios may play a crucial role in deter¬ 
mining the outcome of stellar core collapse in a supernova, 
the fate of a hypermassive, differentially rotating remnant 
of a binary neutron star merger (Baumgarte, Shapiro & 
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Shibata 2000; Shapiro 2000), the generation of gamma-ray 
bursts via the collapse of rotating massive stars to black 
holes (MacFadyen & Woosley 1999), and the collapse of 
supermassive stars to supermassive black holes (see, e.g., 
Rees 1984, Baumgarte & Shapiro 1999, New & Shapiro 
2001 , and references therein.) 

This paper is partitioned as follows: In Section 2 we 
present a brief overview of magnetized Oppenheimer- 
Snyder collapse, and in Section 3 we provide a review of the 
matter and metric fields. In Section 4 we derive analytic 
Newtonian and relativistic solutions for the interior elec¬ 
tromagnetic fields, and in Section 5 we set up Maxwell’s 
equations and boundary conditions for determining the ex¬ 
terior fields. In Section 6 we provide electrodynamic initial 
data, both for the interior and the exterior. In Section 7 we 
present numerical results for the evolution of the exterior 
fields in three different time slices. In particular, in Section 
7.4 we present results in Kerr-Schild coordinates, which al¬ 
low us to follow the evolution to arbitrary late times and 
to track the late-time fall-off of the fields. We briefly sum¬ 
marize our analysis in Section 8 . We also provide two 
Appendices. Appendix A derives a global electromagnetic 
energy conservation criterion, and Appendix B contains 
the finite-difference equations used in our simulations. 

2. MAGNETIZED OPPENHEIMER-SNYDER COLLAPSE: 

OVERVIEW 

In this paper we consider the “simplest” case of mag¬ 
netized stellar collapse to a black hole: the collapse of a 
spherically symmetric, homogeneous dustball, momentar¬ 
ily at rest at t = 0. The matter and gravitational fields 
are described by the exact Oppenheimer-Snyder solution. 

We solve for the electrodynamic fields in a “dynami¬ 
cal Cowling approximation”: the background geometry is 
dynamical and determined by the imploding matter field 
alone; the matter and the geometry (metric) are unaffected 
by the electrodynamic fields, which are evolved in the un¬ 
perturbed background. Our approximation is physically 
applicable whenever the magnetic field is sufficiently weak 
that it satisfies the condition R^/Stt <C GMp/R, where G 
is the gravitational constant, p is the matter density and 
M and R are the mass and radius of the star. Note that 
for a frozen-in B-field, both terms scale the same way with 
radius (i.e. R~'^) during the collapse, so if the inequality 
is satisfied initially, it is satisfied at all later times. 

We seek the relativistic generalization of the solution 
for the electromagnetic field generated inside and outside 
a collapsing, homogeneous, perfectly conducting magne¬ 
tized dust sphere which at < = 0 is threaded by a constant 
interior B-field, matched onto an exterior magnetic dipole 
field. The solution to this problem requires a numerical 
integration of Maxwell’s equations in a dynamical space- 
time in which a black hole forms. While various aspects of 
this problem have been treated previously, including the 
asymptotic field behavior at late times in a static vacuum 
Schwarzschild spacetime, the full dynamical solution has 
never been presented. 

We will divide our time-dependent solution of the elec¬ 
trodynamic fields into two domains: the matter interior 
solution (Section 4) and vacuum exterior solution (Sec¬ 
tion 5). The solutions are joined at the stellar surface by 
applying the junction conditions for an electromagnetic 


field across a moving discontinuity. Hereafter we adopt 
gravitational units whereby G = 1 = c. 


3. THE MATTER AND METRIC FIELDS: THE OS SOLUTION 


The Oppenheimer-Snyder solution (Oppenheimer & 
Snyder 1939) can be derived from the Euler equations, the 
continuity equation and the Einstein field equations. Each 
fluid element follows a radial geodesic in this approxima¬ 
tion. The interior metric is given by the familiar (closed 
Friedmann) line element 

ds^ = —dr"^ + + SIT? (3-1) 


Here r is the (proper) time coordinate, measured from the 
onset of collapse, x is a (Lagrangian) radial coordinate and 
A is related to r through the conformal time parameter 77 , 

A=]^Am(l + cos'q) (3-2) 

T =]^Am{ri + SIT p). (3-3) 


The surface of the star is at some radial coordinate X = Xs 
and the parameter p varies between 0 and tt. The exterior 
metric is given by the Schwarzschild line element, 


ds^ 



dt^ + 



drl + rld^l\ 
(3-4) 


where M is the star’s total mass. The surface of the star is 
at (areal) radius = Rs{t) and follows a radial geodesic 
according to 


Rs = ^i?s(0)(l -k COST?), 

/i?3(o)y/2 ^ 

T = ^ ’ I ip + smp). 

\ 8M J " 


(3-5) 

(3-6) 


Matching the interior and exterior solutions at the surface 
yields 


A — 

-Ci-m — 


sm Xs = 


\ 2M ) 
2M 


(3-7) 


(3-8) 


According to the above equation, Xs must lie in the range 

0 < Xs < t/2. 

The fluid 4-velocity satisfies the geodesic equations. 


= 0 . 


(3-9) 


The rest mass-energy density p measured in the comoving 
frame, remains homogeneous and is given by 


where 


P 

P(0) 


Q""(r), 


Q{t) 

T 


^ = 7:A + cosp) 



1/2 

{p + sin 77 ) 


(3-10) 


(3-11) 

(3-12) 


for 0 < p < TT. The proper time to collapse to a singularity 
at the origin is TcoU = 7 r(i?f ( 0 )/ 8 M)^/^. 
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4. THE INTERIOR ELECTRODYNAMIC FIELDS 


4.2. The Relativistic Solution 


4.1. The Newtonian Solution 


The Newtonian solution describing spherical, homoge¬ 
neous dust collapse has the same functional form as the 
relativistic solution presented in Section 3. In particular, 
the density satisfies 

^ = 0 - 40 , («) 

where Q as a function of Newtonian time t is given by 
(3-11) and (3-12) if r is replaced by t, and the velocity is 


Q{t) 

V = - V 

Q{t) ■ 


(4-2) 


(4-6) 


Here the dot denotes a derivative with respect to time. 

We assume the interior of the collapsing dust ball to be 
perfectly conducting, in which case the electric field van¬ 
ishes in the comoving frame, 

E = 0, (4-3) 

and the magnetic field can be treated in the ideal mag¬ 
netohydrodynamics approximation. It then satisfies the 
induction equation 

atB = Vx(vxB) (4-4) 

and the divergence (or constraint) equation 

V • B = 0. (4-5) 

In index notation in a coordinate basis, equation (4-4) may 
be written 

dtB^ = - v^B*) 

where Di is the covariant derivative associated with the 
spatial metric 7ij, and where 7 is its determinant. In flat 
spacetime with dtj = 0, the last equation can also be ex¬ 
pressed as 

dtB^ = dj{NB^ (4-7) 

where is the magnetic vector density. This 

expression is identical to the relativistic result (4.13) of 
Paper I, which hereafter we will refer to as (1.4.13). In 
terms of the constraint equation (4-5) can be written 

= 0, (4-8) 

which is the same as the relativistic expression (1.4.14). 

In flat-space Cartesian coordinates we have 

* Q i 

and 7 = 1, hence = Hb Equation (4-7) can then be 
expanded as 

dtB^ = v\jB^ + NB\^ - v\jB^ - v^B\j 

= -5^B^ -v^B\. (4-10) 

Q Q 

Since the total derivative following fluid elements is 
d d ^ d 

'(U ~ 


(4-9) 


we find 


with implies 
B%t) 


= -2^H* 
dt Q 


H*(0) 


= = 


(p(o) 


2/3 


RsjO) 

R, 


(4-11) 

(4-12) 

(4-13) 


Here we have used (4-1), and Rg and Rs{0) are the radius 
of the star and its initial value, respectively. 


The relativistic interior solution is most easily evaluated 
in the comoving coordinates of the metric (3-1). In this 
coordinate system the fluid four-velocity is 


or 



(4-14) 


u° = 1, u* = 0. (4-15) 


In these comoving coordinates the relativistic induction 
equation (1.4.13), which is identical to the Newtonian ver¬ 
sion (4-7), reduces to 

Ab* = = 0. (4-16) 

For the metric (3-1) we have 

= A^sin^Xsin^ (4-17) 


where A = A{t) is given by equations (3-11) and (3-12). 
The components of are most easily expressed in an or¬ 
thonormal basis. The orthonormal basis vectors can be 
written in terms of the coordinate basis vectors as 


_ 1 _l d 

1 1 d 

^ A sin X A sin y dO 

1 1 d 

A sin X sin 9 ^ A sin x sin 9 d(j> 


(4-18) 


As a consequence, all orthonormal components are pro¬ 
portional to AH*. Inserting these together with (4-17) into 
(4-16) then yields 


H*A2 = const (4-19) 


or 


H*(0) 




p (Rs{Q)\ 

p(o); \ Rs ) 


(4-20) 


where we have used (3-10). The above equation holds for 
any initial field configuration H*(0). As in the Newtonian 
case, the electrical field H* vanishes in the comoving frame 
under the assumption of perfect conductivity. 


H*(r) =0. 


(4-21) 
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5. THE EXTERIOR ELECTRODYNAMIC FIELDS 


5.1. Maxwell’s equations 

The exterior of the collapsing dust sphere is pure vac¬ 
uum, so that the approximation of ideal magnetohydrody¬ 
namics cannot be applied there. Instead, Maxwell’s equa¬ 
tions have to be solved without approximation. In terms 
of the electric field Ei and the vector potential Ai as mea¬ 
sured by a normal observer n“, Maxwell’s equations can 
be written 

dtA^ = -aEi + P^Aij + P\iAj - d^{a^) (5-1) 

(equation (1.3.22)) and 

dtE^ = 7 - 1/2 (^71/2(7*'7^™ - ,j 

-AiraE + aKE^ + (3 ^E\j - E^ (5-2) 

(equation (1.3.25)). Without loss of generality, we choose 
the source-free electromagnetic Coulomb gauge condition 

$ = 0. (5-3) 

In Section 7 we will present numerical solutions to equa¬ 
tions (5-1) and (5-2) in three different time slicings and 
spatial coordinate systems: Schwarzschild time slicing in 
isotropic coordinates, maximal slicing in isotropic coordi¬ 
nates, and Kerr-Schild time slicing in Kerr-Schild spatial 
coordinates. For the first two spacetime metrics, the met¬ 
ric takes the spherically symmetric isotropic form 

ds^ = —(q;2 — A?l3'^)dt^ 2A^l3drdt + A‘^{dr‘^ + 

(5-4) 

where /? = /?’’, a and A are all functions of r and t alone, 
and K = 0. We will focus on this form here, and will 
generalize to a nonisotropic Kerr-Schild metric in Section 

7.4. 

In this paper we shall consider axisymmetric (dipole) 
field configurations. Complete axisymmetric solutions ex¬ 
ist for which the only non-vanishing components of Ai and 
E* are A^ = A^{t,r,6) and E^ = E^{t,r,0). (Corre¬ 
spondingly, the only non-vanishing components of 5* are 
B^{t,r,9) and B^{t,r,9)). Inserting these quantities to¬ 
gether with (5-3) and (5-4) into equations (5-1) and (5-2) 
yields 

dtA^ = -aE^ + l3drA^ (5-5) 

and 

= ~A^r^ sin9 (( .. + (nr^sin^^^’0 , 0 ) 

+PdrE'^. (5-6) 

The fields and E^ can now be decomposed into mul¬ 
tipoles. A particularly convenient decomposition is 

A^t, r,9) = -J2 r) ^(1 - > (5-7) 


where /i = cosd and where the Pi{iE) are the Legendre 
polynomials of order 1. Inserting this into (5-6) yields 


dtE’^ 


1 

A 3 j .2 



a 






dBi(fx) \ 
dfi ) 


■IddrE'l’, 


where a summation over repeated I is implied, 
satisfy 


A. 

dfj, 



—l{l + 3)Pi, 


equation (5-8) can rewritten as 


dtE'*’ = 


1 (a I \ 

A3^2 ’V, 


1(1 


“,a‘ 


A4 


m 

dfj, 


We now decompose E’l’ as 

E^{t,r,9) = J2 


e\t, r) dPi 
dfi 


The covariant component E^ is then 


dfj, 


Since the 


(5-9) 

-f fddrE^. 
(5-10) 

(5-11) 


(5-12) 


This expression can be inserted into equations (5-5) and 
(5-10), which, after dropping the superscript I in a* and 
e*, yields the coupled set of differential equations in t and 
r for each mode 1: 



a 1(1 + 1) 

A4 j .2 


(5-13) 


dta = aA'^e + f3a^r- (5-14) 

In the following we will specialize to pure dipole fields with 
I = 1. 


5.2. Boundary Conditions 

The dynamical fields a and e satisfy boundary condi¬ 
tions both on the surface of the star and on the outer edge 
of the numerical grid, which we take to be at a large sepa¬ 
ration from the star. On the surface of the star, the fields 
have to be matched to the interior solution (Section 5.2.1), 
except when the stellar interior collapses inside the event 
horizon and is excised from the numerical grid (Section 
7.4). At large distance we impose “outgoing wave bound¬ 
ary conditions” (Section 5.2.2). 

5.2.1. Inner Boundary Conditions 

In a frame comoving with the surface of the star, the 
component of the magnetic field normal on the surface and 
the tangential components of the electric field have to be 
continuous across this surface. To derive boundary condi¬ 
tions for a and e, we therefore first have to boost from our 
computational frame of normal observers n“ to a frame 
comoving with the surface of the star m“. We then have 
to express the transformed junction conditions on the or¬ 
thonormal components of the magnetic and electric fields 
in terms of the expansion coefficients a and e. 

Denote the normal observer frame by E' and the frame 
comoving with the stellar surface by E. Boosting with 
velocity v from E to F' leads to the following relation be¬ 
tween the orthonormal components of the magnetic and 
electric fields in the two frames: 


Efi 

= Ef 

Bf, 

= Bf 

Eg, 

= liEg-vB^) 

Bg, 

= -f(Bg+vE^) 

El, 

4> 

= l{E^ + vBg) 

Bi, 

4> 

= l{B^-vEg) 


(5-15) 
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o(*") _ niout) 
r r 

^(in) _ j^iout) _ Q 


Here 7 is the gamma-factor between the two frames F' 
and F, 

7 = -Uan“ = (1 -w2)-i/2^ (5_16) 

In the comoving frame F, the relevant junction conditions 
for Bf and at the surface are (see, e.g., Eqs. (21.I6Ia-d) 
in Misner, Thorne & Wheeler 1973) 

(5-17) 
(5-18) 

•r •r 

where the electric field component has to vanish because 
of the assumption of MHD in the interior, equation (4-21). 
Combining the Lorentz transformations (5-15) with equa¬ 
tion (5-17) yields the following condition for the magnetic 
field in the frame F' of a normal observer: 

To derive a condition on the electric held , it is useful 

to consider the inverse transformation 

^ _ ^ 4 °“*)). (5_20) 

(p ^ <p' 9' ' ^ ' 

According to (5-18) the left hand side of (5-20) has to van¬ 
ish, which yields the following boundary condition on the 
electric held in F'-. 

(5-21) 

We now express the orthonormal components of the ex¬ 
terior magnetic and electric held in terms of the quantities 
a and e. The magnetic held can be found from 

(5-22) 

(see, e.g., Eq. (1.3.17)). The Levi-Civita tensor can be ex¬ 
pressed as where [ijk] is the completely 

antisymmetric symbol, from which we hnd 

^ (5-23) 

/ aiii 1 / 

and 


B'^ = - 


A^r'^ sin 9 

1 




A. 


(p,r 


(5-24) 


for the metric (5-4). Specializing to pure dipole helds with 
1 = 1, the expansion (5-7) for reduces to 
A(j, = —a sin^ 6, 
and we therefore have 

2 cosd 


B'^ = AB^ =-■ 


and 


B^ = ArB^ = 


e 


A^r 

sin 0 

A^r 


For the exterior electric held E'^ we similarly hnd 
E* = Ar sin 9E’^ = e. 


(5-25) 

(5-26) 

(5-27) 

(5-28) 


Inserting the last three expressions into the boundary 
conditions (5-19) and (5-21) on the helds, we hnally ob¬ 
tain 


a = —A^r‘ 


B. 


(in) 


2 cos 9 


and 


e = , 

A3 ’ 


(5-29) 


(5-30) 


For simulations in Kerr-Schild coordinates, we excise the 
interior of the black holes at late times. In this case, these 
inner boundary conditions have to be replaced with bound¬ 
ary conditions on the surface of the excised region, as we 
will discuss in Section 7.4. 


5.2.2. Outer Boundary Conditions 

In the asymptotically hat region at large distances from 
the star, equations (5-13) and (5-14) combine to give 


0^,tt ^,rr 


l{l + l)a 


(5-31) 


and a similar equation for e. For our adopted I = 1 ini¬ 
tial exterior magnetic dipole held (see Section 6.2), the 
solution to equation (5-31) is characterized by a longitu¬ 
dinal dipole held a ~ 1 /r at early times and an outgoing 
radiation held a ~ f{t — r) at late times. A suitable but 
approximate boundary condition which we adopt to rehect 
this behavior in the wave zone r ^ A ~ 25M is of the form 


e, a 


fit - r) 

r 


(5-32) 


Finite diherence implementations of this equation can be 
found in Appendix B. In fact, we hnd that most of our re¬ 
sults are quite insensitive to the precise form of our outer 
boundary conditions. The reason is that the longitudi¬ 
nal dipole held has a vanishingly small value at our outer 
boundary initially and our integrations terminate before 
any electromagnetic wave ever reaches the outer bound¬ 
ary. 


6. ELECTRODYNAMIC INITIAL DATA 

As initial data for the electromagnetic helds we seek the 
relativistic generalization of the Newtonian data describ¬ 
ing a uniformly magnetized sphere of radius R, matched 
onto an exterior dipole magnetic held. In orthonormal po¬ 
lar coordinates, the interior solution is given by the con¬ 
stant held 

= H(cos0, — sind, 0) (Newtonian interior) (6-1) 

for the held aligned with the z-axis, and the exterior solu¬ 
tion for a pure dipole held is given by 

' BR^ 1 

= —^(cosd, - sind, 0) (Newtonian exterior). 

The junction conditions on the surface of the star require 
that the normal (i.e. radial) component of the held be con¬ 
tinuous across the stellar surface. The tangential discon¬ 
tinuity in implies the presence of a surface current. 


6.1. Relativistic Interior Field 

The relativistic generalization of the Newtonian interior 
solution (6-1) has to satisfy the constraint equation (1.3.9) 

DiB^ = 0 (6-3) 


where Di is the covariant derivative associated with the 
interior metric (3-1). A suitable solution is 

/ \ 2 / \ —2 


= B cos 6 
B^ = -Bsm9 


X 


smx 

X 


smx 


Xs 


sm Xs 
X^ 

sin X 6 


(6-4) 

(6-5) 
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where B is an arbitrary constant in space. This solu¬ 
tion automatically satisfies the magnetic induction equa¬ 
tion during the stellar collapse if B is chosen to vary with 
proper time according to 

/ X 2/3 

(«) 

(see Section 4.2 and equation (4-20)). 

At the surface X = Xs, the radial component which 
can be identified with B'", reduces to 

B^ = BcosO (at surface), (6-7) 

which we will use in the boundary condition (5-29). 

Finally, it is easy to see that the components (6-4) and 
(6-5) reduce to the Newtonian solution (6-1) in the Newto¬ 
nian limit 2M/Rs <C 1. From (3-8), this inequality implies 
that sinxs ^ 1 and hence Xs ~ sinxs- If this holds for 
the surface value Xs, it must also hold for all interior val¬ 
ues X ^ Xs- All the fractions in equations (6-4) and (6-5) 
therefore reduce to unity, so that we recover the Newto¬ 
nian solution (6-1). 


6.2. Relativistic Exterior Field 


An analytic solution for the magnetic field in the exterior 
of a static spherical star as measured by a static observer 
has been given by Wasserman & Shapiro (1983): 


Bf = — 


6/id cos 9 


; (^Xs ln(l 




+ (1 + 


^-1 


( 6 - 8 ) 


Bs = 


Qfid cos 9 


-x^ 

1 — 


(1 - Xs ^)l/27 


(a:,(l-x;i)iAln(l-:r;i) 


(6-9) 


Here is the Schwarzschild (areal) radius, Xs = rg/i^M), 
and Hd is the magnetic dipole moment. Using the junc¬ 
tion condition (5-17), fid can be expressed in terms of the 
interior field by matching the interior and exterior com¬ 
ponents of Bf (equations (6-4) and (6-8)) on the surface. 
This yields 

■ A7V2))-', 

( 6 - 10 ) 

where Rg and Xg are the values of Xg and Xg on the stellar 
surface, and where B is the factor appearing in equation 
(6-4). In the Newtonian limit 1 the solutions (6-8) 

and (6-9) with fid given by (6-10) reduces to the Newto¬ 
nian solution (6-2). 

From equations (6-8) and (6-9) the magnetic potential 
can be found to be 


and therefore, according to equation (5-25), 

(^xl\n{l-x-^) + Xg + ^^ . ( 6 - 12 ) 

We take the above form for the exterior magnetic field 
as initial data in our numerical simulations. We also set 

A* = 0 = e, (6-13) 


which is consistent with our star having zero charge. For 
Schwarzschild slicing and maximal slicing the normal ob¬ 
server is a a static observer at t = 0, a moment of time 
symmetry, so the components of 5* and A* as measured 
by n“ are identical to the values quoted above at the ini¬ 
tial time (see Section 7.4.3 for the case of Kerr-Schild slic¬ 
ing). However, since our dynamical Oppenheimer-Snyder 
spacetime is not the spacetime of a static star, the above 
initial data no longer yields a static magnetic field solu¬ 
tion. However, it does provide a physically plausible field 
for our (arbitrary) initial data, if we assume that star was 
dominated by a dipole field up to the moment of collapse. 


7. EXTERIOR ELECTROMAGNETIC FIELD SOLUTION IN 
DIFFERENT TIME SLICINGS 

We now integrate the coupled equations (5-13) and (5- 
14) for the exterior electric field and the magnetic po¬ 
tential A(f,, given the interior solution and initial data as 
in Sections 4 and 6. 

This problem provides an ideal laboratory for exper¬ 
imenting with various coordinate gauge choices for per¬ 
forming numerical integrations in strong gravitational 
fields. We will compare three different time slicings 
for treating the exterior evolution, namely Schwarzschild, 
maximal and Kerr-Schild. Historically, the earliest nu¬ 
merical studies traditionally worked in the Schwarzschild 
gauge when dealing with spherical, asymptotically flat 
spacetimes. But in this gauge the time coordinate ap¬ 
proaches infinity as the stellar surface reaches the black 
hole event horizon at areal radius Vg = 2M; light cones 
pinch off near the horizon and the resulting coordinate 
singularities make integration of the electromagnetic field 
equations to arbitrarily late times difficult in this gauge. 
Maximal time slicing has often been adopted in recent nu¬ 
merical work as a “singularity avoiding” gauge; spatial 
slices cover the entire matter interior as well as the vacuum 
exterior and penetrate the horizon. But “grid stretching” 
along the black hole throat (see, e.g., Shapiro & Teukolsky 
1986) as the stellar surface approaches a limiting surface 
inside the horizon makes the solution increasingly inac¬ 
curate at late times. Maximal time slicing requires the 
numerical solution of a second order elliptic equation for 
the lapse. While boundary conditions on this equation are 
straightforward to impose at the stellar origin and at large 
distances from the matter field, they are not at all simple 
to impose at a black hole horizon. Consequently, it is not 
straightforward in this gauge to “excise” the black hole in¬ 
terior from the computational grid, even though this region 
(where inaccuracies are growing) is causally disconnected 
from the exterior. These problems can be avoided in Kerr- 
Schild slicing, where the lapse is nonzero on the horizon 
(unlike Schwarzschild but like maximal) and analytic (like 
Schwarzschild but unlike maximal). These features allow 
for horizon penetration and “black hole excision” (Unruh, 
unpublished; Thornburg 1987; Seidel & Suen 1992) to fol¬ 
low the evolution of the electromagnetic fields outside the 
black hole to arbitrary late times. 
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7.1. Analytic Precursor: Quasi-Static Approach 


In the earliest attempts to calculate the exterior mag¬ 
netic field of a magnetized star undergoing stellar collapse 
it was assumed that the exterior field is a (longitudinal) 
dipole field which changes quasi-statically during the entire 
collapse. In this approximation, the electric field and the 
time derivative of the magnetic field is neglected, and the 
exterior magnetic field is given, at all times, by equations 
(6-8) - (6-9). The magnetic field strength is determined 
by the magnetic dipole moment fj,d, which can be found 
from equation (6-10). As the stellar surface Rs approaches 
the horizon, the magnetic dipole moment approaches zero 
according to 


1 


/i oc 


ln(l - 2M/Rs) 


as R, 


2M. 


(7-1) 


During this late phase of the collapse we have 


1 - 


2M 


(X exp(—t/4M) as Rs 


2M 


(7-2) 


(see, e.g. Problem 15.11 in Lightman et al. 1975). To¬ 
gether, these imply 

pL oc t~^ as Rs 2M, (7-3) 


suggesting that the exterior magnetic field should decay 
with t~^ at late times (see, e.g., Ginzburg & Ozernoy 1965; 
Anderson & Cohen 1970; Zeldovich & Novikov 1971). This 
result is incorrect, since the time-changing magnetic dipole 
in the collapsing star generates a transverse electromag¬ 
netic field {E and B), which dominates the late-time be¬ 
havior but is ignored in this quasi-static approximation. 

The correct decay rate at late times of an initially static 
dipole electromagnetic radiation field outside a black hole 
is actually ^-(2^+2) (Price 1972a,b; also pointed out by 
Thorne 1971). Our numerical integrations below confirm 
this result for / = 1. In addition, we are able to obtain 
the complete interior and exterior numerical solution over 
the entire collapse of a magnetic star to a black hole and 
show the transition from a quasi-static longitudinal to a 
dynamic transverse radiation field. 


7.2. Schwarzschild Slicing 

In Schwarzschild time slicing and isotropic coordinates, 
the exterior metric takes the form (5-4) with the lapse a 
given by 

, 1 - M/W 

1 + M/(2r) ’ ' ' 

the conformal factor A given by 

+ (7.5) 

and with the shift /3 vanishing identically, 

P = 0. (7-6) 

In the above, the isotropic radius r is related to the areal 
radius rs by 

mV 


rs=r\l + -) 


(7-7) 


and 


i (r. - M + {rsirs - 2M)f/^) . (7-8) 


The above quantities appear in the evolution equations 
(5-13) for e and (5-14) a. 

Analytic expressions for the areal radius Rs of the stel¬ 
lar surface as a function of proper time t have been given 
in Section 3 in terms of the parameter ry. The isotropic 
radius can be found by substituting the areal radius into 
(7-8). The Schwarzschild coordinate time t can also be 
found in terms of ry according to 


t = 2M\n 


+2M 


f (R,(0)/(2M) - 1)1/2 ^ tan(ry/2) \ 
V(i?.(0)/(2M)-1)1/2-tan(,y/2)j 


(Rsio) 

V 2M J 



Rs{0) 

AM 


{rj -\- sin?y) 


7-9) 


(see, e.g. Misner, Thorne & Wheeler 1973, eqn. 32.10b). 
Note that we have t/M ^ 00 as the stellar surface ap¬ 
proaches the event horizon, Rs 2M. 

The velocity v of the surface with respect to a normal 
observer is 


f2M/Rs-2M/RsiO)y^^ 
V 1 - 2M/Rs{0) J 


(7-10) 


This quantity enters the inner boundary condition for e, 
equation (5-30). 

We have solved equations (5-13) and (5-14) in 
Schwarzschild slicing with an implicit finite differencing 
scheme (see Appendix B). We regrid after each timestep 
so that the stellar surface is always located half-way be¬ 
tween the two first grid-points (which makes the imple¬ 
mentation of the boundary conditions (5-29) and (5-30) 
particularly simple). The advantage of implicit schemes is 
that no Courant condition has to be satisfied to maintain 
stability of the numerical integration. Because of the lin¬ 
earity of Maxwell’s equations the initial field strength Bq 
is arbitrary (provided it is dynamically unimportant) and 
can, without loss of generality, be chosen to be unity. 

In Figures 1 and 2 we show numerical results for a star 
collapsing from rest from an initial radius of 77^(0) = 4M. 
In these snapshots we show the collapsing star together 
with exterior magnetic field lines. In axisymmetry, these 
field lines coincide with contours of constant A^ (Zhang 
1989), which enables us to draw them easily. At t = 0 
all field lines are attached to fluid interior and emerge 
from the stellar surface. At later times, as the surface 
approaches the event horizon, field lines that emerge from 
the surface are crushed against the star and become in¬ 
creasingly tangential (Anderson & Cohen 1970). At larger 
separations, field lines detach from the star during the 
collapse. In the radiation zone, these disjoint lines propa¬ 
gate outwards with the speed of light. Thus, early on, we 
find that the exterior longitudinal B field evolves quasi- 
statically to match the growth of the frozen-in interior field 
(see equation (4-20)), while the exterior E field remains 
small. At late times, the exterior longitudinal fields are 
radiated away and the dominant B and E fields become 
outgoing, decaying, transverse electromagnetic waves. 

To calibrate the accuracy of the numerical integration 
we check the global conservation of magnetic energy as 
outlined in Appendix A. We show results of this energy 
conservation test in Figure 3. Clearly, the accuracy dete¬ 
riorates soon after the surface of the star approaches the 
event horizon at around t ~ 20M. This effect is not sur¬ 
prising. The vanishing of the lapse (7-4) on the horizon 
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t/M = 0.00 


t/M = 8.9 




t/M = 13.3 


t/M = 13.4 




t/M = 15.9 t/M = 17.7 



Fig. 1.— Snapshots of the exterior dipole magnetic field lines at select Schwarzschild time slices for a star which collapses from rest at 
an initial areal radius Rs = AM. Points are plotted in a meridional plane using areal radii. The light shaded sphere covers the matter 
interior. The initial growth of the longitudinal field due to flux freezing in the interior is ultimately followed by a outward burst of transverse 
electromagnetic radiation as the star approaches the black horizon at Vs = 2M. Soon after the last snapshot at t = 17.7M the quality of the 
numerical integration deteriorates in this gauge. See text for details. 


causes the local light cone to “pinch off” and the (coordi¬ 
nate) speed of light to go to zero. As a consequence, incom¬ 
ing radiation piles up in front of the horizon and develops 
smaller and smaller spatial structures, which ultimately 
become smaller than any fixed grid resolution (cf. Rezzolla 
et at. 1998). 

It is evident that Schwarzschild coordinates are not suit¬ 
able for calculating the late-time radiation as measured by 
a distant observer, and, hence, for determining the late¬ 


time electromagnetic field decay rate. This observation 
motivates using a coordinate system that is regular on the 
event horizon and extends smoothly into the black hole 
interior. 

7.3. Maximal Slicing: Interior and Exterior Solution 

We now follow the same collapse but adopt maximal 
time slicing and isotropic (minimal shear) spatial coordi¬ 
nates. This coordinate system extends regularly into the 
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t/M = 0.00 


t/M = 8.9 




t/M = 13.3 


t/M = 13.4 




t/M = 15.9 


t/M = 17.7 




Fig. 2.— Blow-up of the same collapse depicted in Figure 1. Note how in this gauge the exterior field lines which are linked to the interior 
are crushed tangentially toward the surface as it approaches the horizon. Field lines which pinch off and become disjoint from the interior 
propagate outward as a dipole electromagnetic wave. 


black hole. Maximal slicing is imposed by requiring the 
trace of the extrinsic curvature to vanish at all times: 


K = 0 = 


dK 

dt 


(7-11) 


Taking the trace of equation (1.2.8) then yields an elliptic 
equation for a, which in general has to be solved numer¬ 
ically. Maximal slices of Oppenheimer-Snyder spacetime 
in isotropic coordinates may be determined by solving or¬ 
dinary differential equations. Solutions are presented in 
Petrich, Shapiro & Teukolsky (1985) and provide a, /3 and 


A, which we substitute in equations (5-13) and (5-14) for 
the evolution of e and a. These coordinates cover both the 
interior and the exterior of the star, so that we can map 
out the interior as well as the exterior field on these slices. 

The Lorentz factor between a normal observer and a co¬ 
moving observer at the stellar surface 7 = —Uan°' can be 
computed from the numerical solution presented in Pet- 
rich, Shapiro & Teukolsky (1985) (i.e., 7 = 1/(1 —T^) and 
V = |r|, where T is given by equation (6.11) of that pa¬ 
per). This is the velocity that is needed in the boundary 
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Fig. 3.— Electromagnetic energy conservation in Schwarzschild time slicing. The outer boundary is located at i^out = lOOM, the exterior 
of the star is covered with 5000 gridpoints, and the Courant factor is unity. The solid line denotes the total energy E, which should be 
conserved (see equation (A8)). The dotted line is the energy S{t) contained between the radii r\ = 1.01i?s(0) = 4.04M and r 2 = SOM; the 
short dashed line is the integrated flux across ri, J7(ri), and the long dashed line is the integrated flux across r 2 , All energies are 

normalized to the initial energy Eq. Note that the quality of the numerical integration starts to deteriorate at about t 20M, when the 
stellar surface and the ingoing electromagnetic waves approach the event horizon at Vs = 2M. 


condition (5-30). 

We numerically solve equations (5-13) and (5-14) on 
maximal slices using implicit finite difference equations 
that are identical in form to those used on Schwarzschild 
slices (see Appendix B). 

In Figures 4 and 5 we show snapshots of the collapse in 
maximal slicing. Unlike Schwarzschild slices and the Kerr- 
Schild slices below, maximal slicing also covers the matter 
interior, which allows us to trace the interior magnetic 
field lines as well as locate the apparent and event hori¬ 
zons in these plots. We can also identify the limit surface 
at Rs « I.5M which the matter surface approaches at late 
times (Petrich, Shapiro & Teukolsky 1985). Note that in 
the interior, a field line must connect the same Lagrangian 
fluid element at all times to satisfy the “frozen-in” require¬ 
ment of MHD. This criterion provides a check on our field 
line tracer. For comoving observers in the interior, the 
electric field vanishes, but as viewed by a normal observer 
in maximal slicing the electric field is non-vanishing in 
general (see Section 5.2.1 for the transformation). 

In maximal slicing the metric is not static, so that we 
cannot simply apply the electromagnetic energy conserva¬ 
tion check developed in Appendix A. 

Maximal slicing has some very desirable properties, but 
at late times grid stretching effects make the computa¬ 
tion even of the background models increasingly difficult 
(Shapiro & Teukolsky 1986). We terminated the calcula¬ 
tion shortly after t = 20M. 

7.4. Kerr-Schild Slicing: Black Hole Excision 

The problems of maximal slicing can be avoided by us¬ 
ing Kerr-Schild (or ingoing Eddington-Finkelstein) coordi¬ 
nates, which are analytic but smoothly extend across the 
horizon. This slicing is not “singularity avoiding”, so that 
it requires that the interior of a black hole be excised from 


the numerical grid. Adopting horizon-penetrating coordi¬ 
nate systems together with black hole excision is currently 
the most promising approach to dynamical simulations of 
black holes. In this calculation we apply this approach 
to the collapse of a magnetized star. Since the metric 
in Kerr-Schild coordinates (equation (7-12) below) takes 
a form that is different from (5-4), several parts of the 
above calculations, including the form of Maxwell’s equa¬ 
tions, boundary conditions and initial data, have to be 
modified. 


7.4.1. Kerr-Schild coordinates 

In Kerr-Schild coordinates, the metric takes the form 

,2 2M, , 2 4M , , 2M, , 2 2,^2 

ds^ = —(1- )dt^ H- dtdrs + (IH- )dr^g r^dfl^, 

rs rs rs 

(7-12) 

where is the Schwarzschild (or areal) radius. From the 
metric we can identify the lapse as 

a= f - ' , (7-13) 

\rsE2M) ' 

the radial shift as 

2 M 

(3 = -—, (7-14 

rs-\-2M ^ ' 

and the diagonal components of the spatial metric as 

2M 

= 1 + - (7-15) 


169 — . 2 ~ '^s 

sm 0 


(7-16) 


(see, e.g., Lehner et al. 2000). The extrinsic curvature can 
then be computed from equation (1.2.9), which yields 


2M, 

= -^(rs -I- M)a 


(7-17) 
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t/M = 17.0 
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Fig. 4.— Snapshots of the interior and exterior magnetic field lines at select maximal time slices for the same collapse depicted in Figure 1. 
The star again collapses from rest at an initial areal radius /?s(0) = 4M. Points are plotted in a meridional plane using areal radii. The light 
shaded sphere covers the matter interior; the black shaded sphere covers the region inside the event horizon. Soon after the last snapshot at 
t = 20.7 grid stretching causes the quality of the numerical integration to deteriorate in this gauge. See text for details. 


and 

Kee = = 2Ma. (7-18) 

sin 0 

The trace of the extrinsic curvature is 

2M 1 -f 3M/rs _2Ma rs + 3M 
r2 (1-h 2M/r,)3/2 “ r2 r, + 2M' 

(7-19) 


The surface of the collapsing star follows a radial 
geodesic, which can be determined as follows. The areal 
radius as a function of proper time has to satisfy the same 
equations, 

Rs = (1 + cos? 7 ) 



7.4.2. Radial geodesics in Kerr-Schild coordinates 


8M 


{rj + sinry) 


(7-20) 
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t/M = 0.00 
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t/M = 11.9 
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t/M = 17.0 


t/M = 20.7 




Fig. 5.— Blow-up of the same collapse depicted in Figure 4. Note how in this gauge the stellar surface approaches a limit surface at 
Vs ^ 1.5M at late times. The horizon grows to its final value Vs = 2M once all the matter crosses inside this radius. 


as before (see Section 3), since both quantities are 
gauge invariant parameters. Here i?s(0) is the initial 
Schwarzschild radius at which the surface is momentar¬ 
ily at rest. An equation for the Kerr-Schild coordinate 
time t can be derived from the conservation of pt 


Combining this equation with the normalization condition 
for the four-velocity, UaU°' = — 1, yields 


dvs 

dr 



(7-23) 


Pt = gttp^ + gtr,p''‘ = const = -E, 


or 


(7-21) 


Evaluating this result for the surface at t = 0, when 
dRs/dr = 0, determines the energy 


2M\ dt 
Ts ) dr 


2M dvs _ ~ 
Ts dr 


(7-22) 


E = 



2M \ 

'iuo)) 


1/2 


(7-24) 
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Fig. 6.— Electromagnetic energy conservation in Kerr-Schild time slicing. The outer boundary is located at Rout = lOOM, the exterior of 
the star is covered with 5000 gridpoints, and the Courant factor is 0.5. The solid line denotes the total energy E, which should be conserved 
(see equation (A8)). The dotted line is the energy S{t) contained between the radii ri = 1.01Rs{0) = 4.04M and r 2 = SOM; the short dashed 
line is the integrated flux across ri, J7(ri), and the long dashed line is the integrated flux across r 2 , Note that approximately 80% of 

the energy enters the black hole, while only 20% is radiated away to inflnity. 


Inserting the last two equations into (7-22) then yields 
dt (1 - 2M/i?,(0))i/2 2M (2M/R,-2M/R,(0))1/2 
dT~ 1 - 2M/Rs ~ Rs 1 - 2M/Rs 

(7-25) 

Interestingly, this equation can be integrated analytically 
to yield 


, = 2m(AA_i 

V 2M 


-{ij + sinr])] (7-26) 


-t 2 Mln 


1 -L _ 

^ ^ 2M 


1(3^ _ 

V2M 


The coordinate time tf for infall from = i?s(0) > 2M 
to the singularity at rg = 0 is therefore 


G = 27rM 1 - 5 ^^ - 2 Mln 

^ \ AM ) \ 2M J \ 2M 

(7-27) 

This time is evidently finite, which reflects the fact that 
Kerr-Schild coordinates are horizon penetrating. This also 
demonstrates the necessity of excising the black hole inte¬ 
rior, since otherwise the central singularity would be en¬ 
countered after t = tf. 

The Lorentz factor corresponding to a boost from an ob¬ 
server comoving with the stellar surface to a local normal 
observer can be found from 

7 = -naU°' = au* = ( 1- — ] . (7-28) 


The relative velocity v of the normal observer as measured 
by a comoving observer at the surface is then given by 

v = -^. (7-29) 

iXs 

A minus sign appears in equation (7-29) because the shift 
drives the radial velocity of normal observer inward more 
rapidly than the surface. 


7.4.3. Exterior initial data in Kerr-Schild coordinates 

Given the boost parameters 7 and n, we can transform 
the magnetic dipole initial data of Section 6.2 into Kerr- 
Schild coordinates. Applying the transformation rules (5- 
15) we find 

Dr _ Dr 

— -^WS 

4s = 74s 

4s = ^ 

4s = 4s = 0 

4s “ 7(^ws ~ ^-®ws) 

4s = 7(4s+^4s) 

"Here the subscript KS refers to Kerr-Schild data, and WS 
to the Wasserman-Shapiro exterior magnetic dipole solu¬ 
tion presented in Section 6.2. We note that even at t = 0 
the electric field is non-vanishing as seen by a normal ob¬ 
server in Kerr-Schild coordinates, since a normal observer 
is moving with respect to a static observer. 


7.4.4. Maxwell’s equations in Kerr-Schild coordinates 

In Section 5.1 we evaluated Maxwell’s equations (5-1) 
and (5-2) for a metric of the form (5-4). Since the Kerr- 
Schild metric (7-12) has a slightly different form, we have 
to re-evaluate Maxwell’s equations. 

Our derivation closely follows that of Section 5.1. We 
insert the metric coefficients of Section 7.4.1 into equa¬ 
tions (5-1) and (5-2), choose the Coulomb gauge condition 
by setting <I> = 0 , ensure axisymmetric fields by assuming 
that Af, = A^{t,r,6) and = E‘l’{t,r,6) are the only 


(7-30) 

= 0 

= 7^'4s■ 
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Fig. 7.— Snapshots of the exterior magnetic field lines at select Kerr-Schild time slices for the same collapse depicted in Figure 1. The star 
again collapses from rest at an initial areal radius Rs{0) = 4M. Points are plotted in a meridional plane using areal radii. The white shaded 
sphere covers the matter interior; the black shaded area covers the region inside the event horizon; the grey shaded area covers the region 
inside rs = M excised from the numerical grid once the matter passes inside. In this gauge, using excision, we are able to integrate reliably 
to arbitrary late times. See text for details. 


non-vanishing components, and finally use the expansions 
(5-7) and (5-11). The resulting equations for e and a are 

o ^ 2 \ l{l 1) rt 2 f ^ \ 

Ote = a (a a,rj ^ ^ - 2 — ^ [ ~ I 

’ ^ \^s / ,rg 

(7-31) 

and 

dta = ae + , (7-32) 

where the metric coefficients a and (3 and the trace of the 


extrinsic curvature K are given in Section 7.4.1. 

In terms of e and a, the components of the electric and 
magnetic fields are given by 


E* 


sin0 
-e 

Ta 

2 cos 9 


(7-33) 


(7-34) 
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Fig. 8.— Blow-up of the same collapse depicted in Figure 7. By the end of the integration, all exterior electromagnetic fields in the vicinity 
of the black hole have been captured or radiated away. 


= 


siii0 




(7-35) 


7.4.5. Boundary Conditions 


For the early part of the evolution, the inner boundary 
conditions result from junction conditions on the electro¬ 
magnetic fields at the stellar surface. The derivation of 
these conditions follows that of Section 5.2.1 and results 
in 


tBi 


(7-36) 


and 


V 

® ^ (1 -t 2Mlrsy/'^ 


(7-37) 


Once the stellar surface has passed through some radius 
i?ex < 2M inside the horizon, we switch from the junction 
condition on the stellar surface to an excision boundary 
condition at Rex- We typically set Rex = M, but the 
results are insensitive to the choice provided Rex is well 
inside the horizon at 2M. Once we have switched to an 
excision boundary condition, we fix the numerical grid so 


a = — 


2 cosd 
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Fig. 9.— Far-zone view of the same same collapse depicted in Figures 7 and 8, but now extending to very late times. Note the transformation 
of the dipole magnetic field from a quasi-static longitudinal to a transverse electromagnetic wave. 


that Rex is located half-way between the first two grid- 
points. 

We have experimented with a number of different ex¬ 
cision boundary conditions, but found, as expected, that 
the exact implementation had hardly any effect on the re¬ 
sults in the exterior of the black hole, as long as Rex is 
far enough inside the event horizon. Following a similar 
approach by Alcubierre & Briigmann (2001) for dynamical 
black hole simulations, we settled on simply copying the 
values of e and a from the second gridpoint to the first at 


every timestep, which is equivalent to imposing 

e,rs = 0,,r, = 0 (7-38) 

(see equations (B13) and (B14)). 

7.4.6. Numerical results 

We have implemented the evolution equations (7-31) 
and (7-32) using both implicit and explicit finite difference 
methods (see Appendix B for the finite difference equa¬ 
tions). The results presented here were obtained with the 
explicit scheme. 
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Fig. 10.— The absolute value of the radial magnetic field |B’'| at r = 20M as a function of time. For this simulation, the outer boundary 
was placed at 2000M, the numerical grid consisted of 6000 radial gridpoints, and the Courant factor was 0.5. The dashed line is a power law 
which approximates the late-time fall-off extremely well. 




Fig. 11.— Same as Fig. 10, except for \B^\. 


Since the Kerr-Schild metric (7-12) is independent of 
time, we can calibrate the performance of the code with 
the same global electromagnetic energy conservation test 
as we did in Schwarzschild coordinates (Section 7.2). The 
necessary expressions can be found in Appendix A. We 
show results for the energy conservation check in Figure 6. 
It is quite obvious that this integration works much better 
than that in Schwarzschild coordinates (cf. Figure 3), due 
to regularity across the horizon. We can follow the elec¬ 
tromagnetic fields to arbitrary late times, and can track 
the late-time behavior reliably. 

Figure 6 also shows that about 80% of the energy con¬ 
tained in the initial magnetic field is absorbed by the black 
hole, while about 20% of the radiation is emitted in an 


electromagnetic wavepacket. 

In Figures 7 through 9 we show snapshots of magnetic 
field lines for a star collapsing from Rg (0) = 4M as viewed 
in three different regions. It can be seen very clearly how 
the collapse of the magnetized star emits a pulse of electro¬ 
magnetic radiation, which travels to infinity at the speed 
of light. Similar behavior is seen in plots of the electric 
field at late times. 

We can now evaluate the late-time fall-off of the elec¬ 
tromagnetic fields. In Figures 10 through 12 we show the 
values of \B^\ and \E‘^\ at r = 20M as a function 
of time. As expected, the fields are nearly constant up to 
t = 2QM. Once the collapse gets underway, the exterior 
longitudinal field starts to evolve to match the increas- 
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Fig. 12.— Same as Fig. 10, except for \E^\. 


ing interior frozen-in field. Soon after t = 20M, the first 
electromagnetic pulse induced by the collapse passes the 
observer at r = 20M. Sign changes in the field appear as 
down-ward spikes in the plots. After t ~ lOOM, the exte¬ 
rior fields are predominantly characterized by a decaying 
dipole with a power-law fall-off with time. In Figures 10 
through 12 we have plotted the power-laws curves vary¬ 
ing like to which the fields approach for late times. 
This finding for the exterior fields agrees with the general 
results of Price (1972a,b; cf. Thorne, 1971) who showed 
that multipole electromagnetic fields scattering off a static 
Schwarzschild black hole should decay as 

Measurement of the wavelength of the late-time outgo¬ 
ing electromagnetic waves provides further verification of 
our numerical integrations. We find that the wavelength 
satisfies A « 24.7M. This is the value of a quasi-normal 
mode for electromagnetic dipole radiation scattering off 
a static Schwarzschild black hole (Ferrari & Mashhoon 
1984). This finding shows that the late-time radiation is 
ringing radiation that has a frequency determined by the 
mass of the hole and is independent of the initial electro¬ 
magnetic field profile. 

8. SUMMARY 

We have solved for the evolution in general relativity of 
the electromagnetic field of a magnetic star that collapses 
from rest to a Schwarzschild black hole. We adopted a “dy¬ 
namical Cowling” model by which the matter and gravi¬ 
tational field are prescribed by the Oppenheimer-Snyder 
solution for spherical collapse. Starting from Maxwell’s 
equations in 3 -I- 1 form, we solved for the magnetic and 
electric fields both in the vacuum exterior and the stellar 
interior. We assumed that the initial star was threaded by 
a dipole magnetic field and that the interior behaves as an 
MHD fluid. We showed how the electromagnetic junction 
conditions could be applied at the stellar surface in order 
that the fields in the MHD interior, which were determined 
analytically, could be matched to the fields in the vacuum 


exterior. The exterior fields were evolved numerically. 

To gain experience with solving Maxwell’s equations nu¬ 
merically in a dynamical spacetime containing both mat¬ 
ter and a black hole, we solved the above problem in three 
different coordinate gauges: Schwarzschild, maximal and 
Kerr-Schild. While the first two coordinate systems have 
singularity avoiding properties, they both lead to growing 
numerical inaccuracies in the exterior as the surface ap¬ 
proaches the horizon (in the case of Schwarzschild time 
slicing) or a limit surface inside the horizon (in the case of 
maximal slicing) at late times. Kerr-Schild slicing is hori¬ 
zon penetrating and does not suffer from grid stretching; 
accordingly this slicing enables us to apply excision bound¬ 
ary conditions once the star collapses inside the horizon. 
Adopting Kerr-Schild slicing coordinates with black hole 
excision allowed us to integrate Maxwell’s equations to 
arbitrary late times. We followed the transition in the ex¬ 
terior of a quasi-static, longitudinal magnetic field, which 
evolved during the collapse in to match the growing frozen- 
in interior field, to a transverse electromagnetic wave. By 
the end of the simulation, the electromagnetic field energy 
in the near zone outside the black hole had either been 
captured by the hole or radiated away to large distances, 
in accord with the “no-hair theorem”. At late times we re¬ 
covered the expected power-law decay of the dipole fields 
with time, as well as the expected wavelength of the ring- 
down radiation. 

There exist a number of very important astrophysi- 
cal problems which require the solution of the coupled 
Maxwell-Einstein-MHD equations in a strong, dynamical 
gravitational field. For example, learning the outcome 
of rotating stellar core or supermassive star collapse, the 
mechanism for gamma-ray bursts, and/or the fate of the 
remnant of a binary neutron star merger may all require 
solving this coupled system of equations. The solution 
we have presented for our simple collapse scenario should 
serve as a preliminary guide to the design and testing of 
more sophisticated codes capable of handling these more 
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complicated scenarios. 
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APPENDIX 

GLOBAL ELECTROMAGNETIC ENERGY CONSERVATION 
General Discussion 

The conservation of energy can be verified very easily if the exterior metric is independent of time, so that 


is a Killing vector. In that case, 
is a conserved current and satisfies 


^ dt 


= 0 , 


dV 


(Al) 

(A2) 

(A3) 


where dV denotes the 3-surface enclosing the 4-volume V. If the volume V is confined by the radii ri and r 2 and the 
times ti and t 2 , the conservation law can be written as 

/•r’2 /•TT /*27r d2 


/ r=ri J 9=0 J 0=0 

/•i2 / 


J^^/^drd0d(j) 


^277 


Jt—ti J J(p—O 

Defining the energy contained between ri and r 2 at time t as 


rV^dtdOdcj) 


fr2 /*7r pZTT 

E{t)= / / J^^/^drdOdcj) 

Jr—ri J9—0 J4>—0 

and the integrated flux across the radius r as 

pt2 PTT ^277 

J(r) = / / / r^/^dtdedij) 

Jt=ti 9 9=0 70=0 


the conservation law can be rewritten as 


(A4) 

(A5) 

(A6) 

S{t2)-Siti)=J{r2)-J{ri). (A7) 

The meaning of this statement is quite evident: the difference in the energy contained between ri and r 2 at two different 
times has to equal the net energy flux entering the region, integrated over the time interval. An alternative way of writing 
the conservation law is 

E/Eo = {£{t2) - J{r2) + JirO) IEq = 1, (A8) 

where Eq = E(ti) is the initial energy between ri and r 2 - We use equation (A8) in plotting our results in the text. 

Since the exterior metric is independent of time in both Schwarzschild coordinates (Section 7.2) and Kerr-Schild coor¬ 
dinates (Section 7.4), energy conservation can be verified quite easily in both of those coordinate systems. To do so, the 
energy (A5) and the flux (A6) have to be expressed in terms of the relevant metric quantities, which we do separately in 
the two following Sections. A similar analysis was used in Appendix A of Pavlidou, Tassis, Baumgarte & Shapiro (2000) 
for the energy conservation of a scalar field in the static spacetime of a spherical, equilibrium star. 

Schwarzschild Slicing 

In isotropic Schwarzschild coordinates the metric takes the form (5-4) with a, A and f3 given by (7-4), (7-5) and (7-6). 
The determinant of the metric is 

= aA^r"^ sin0, (A9) 

and the fluxes and J'’ in (A5) and (A6) can be expressed as 


jO rpO rpO rpOO ^ ^T7'2 

J =To=Tg = -T =-^(E 
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Here we have used the results of Section 7 of Paper I to evaluate the source terms. The fields E and B can now be 
expressed in terms of e and a using equations (5-26), (5-27) and (5-28) for an I = 1 dipole field. The angular integrals in 
equation (A5) and (A6) can be carried out analytically, yielding 

N 1 P /^2 2 4 a2 2 (a^)2' 


A3 


(A12) 


and 


J{t) = -a^A / dt{a^re). 
3 Jti 


(A13) 


These expressions can now be inserted into equation (A8) to verify energy conservation in Schwarzschild coordinates. 


Kerr-Schild coordinates 

In Kerr-Schild coordinates, the determinant of the metric is 

= rl sin 9. 


(A14) 


The fluxes J° and in (A5) and (A6) are easiest evaluated by using the form (1.7.8) for the electromagnetic stress-energy 
tensor We then Hnd 

(A15) 


4^J° = 4^r0 = ln°no{E,E^ + B,B^) + n°eoi^E^B'^. 


Here eo/m = —fieir^m (see (1.3.19)) and restricting the summations to the only non-vanishing components E^, B'' and B^ 
yields 

47rJ° = -i(E^ + B^) - 1 + E^B^. (A16) 

2^ ’ l + rj2M ^ ’ 


■rsl2M 

The radial component can be evaluated in a similar way, which leads to 


47rJ'’ = 
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l + r,/(2M) 

1 -h 2M/r. 


({E^f + (H®)' 
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e^bK 


,(l + r,/(2M))2 l + 2M/r,; 

Specializing to / = 1 dipole fields, the angular integrals in equation (A5) and (A6) can again be carried out analytically, 
which yields 

r ’’2 _ /9 ^=2 A „2 
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and 
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1 




.(l + r,/(2M))2 (l + 2M/r,)3/2, 

These expressions can be inserted into equation (A8) to verify energy conservation in Kerr-Schild coordinates. 


FINITE-DIFFERENCE EQUATIONS 

To finite difference the coupled evolution equations (5-13) and (5-14) for Schwarzschild or maximal slicing (or (7-31) 
and (7-32) in the Kerr-Schild case) we introduce a uniform radial grid = 1,2,...,/. We regrid the radial gridpoints 
after each timestep, using linear interpolation between the old gridpoints, so that the stellar surface Rs is always half-way 
between the first two gridpoints at (ri -|- r2)/2. In the Kerr-Schild simulations this is no longer necessary once the interior 
of the black hole has been excised and the inner boundary condition is imposed at a constant excision radius i?ex. All 
variables are defined at the gridpoint r^. Note that our numerical grid never extends to r = 0. In the following we use 
At = - r and Ar = n+i - r,. 


Schwarzschild Coordinates and Maximal Slicing 

For isotropic Schwarzschild coordinates (Section 7.2) and maximal slicing (Section 7.3) we adopt an implicit evolution 
scheme. It is convenient to use the abbreviations 

B = 5 (B1) 
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C = A^. 

Finite difference representations of equations (5-13) and (5-14) can then be written as 


- a” 


_ j^n+lQin+l^n+l 




e”+^ - e? 


j DK+l („n+l 

1+1/2 (®i-|-l ®i-l 


/(H-jO^n+1 


Here the midpoint values are computed as linear interpolations 

differenced the shift term implicitly in (B3), but explicitly in (B4). 

The boundary condition (5-29) is implemented as 


^)/2. Note also that we have finite- 


^ =-(Ay,V3/.)^|, (B5) 

where the interior field strength B is given by equation (6-6). The outer boundary condition (5-32) is finite differenced as 

At Ar n ^ ^ 

Substituting the finite difference equation (B4) into (B3) and using equations (B5) and (B6) at the inner and outer 
boundaries yields a single tridiagonal matrix equation for For this approach to work we had to finite-difference the 

shift terms in (B4) explicitly, because these would have otherwise introduced terms into the equation for We 

invert the tridiagonal matrix at each time step to obtain the new values New values can then be found by 

using equation (B4) again together with the inner boundary condition (5-30) 

e"+i+e?+i t’s/V a?+^-a?+^ 


and an outer boundary condition equivalent to (B6). 

To choose the size for successive timesteps we set 

At = C min (Ar), (B8) 

where we have introduced a Courant factor C in the above equation. Implicit differencing imposes no restrictions on C, 
but accuracy usually requires the choice C < 1. 

Kerr-Schild Slicing 

The evolution equations (7-31) and (7-32) in Kerr-Schild slicing can be finite differenced in a very similar fashion to 
those in Schwarzschild coordinates or maximal slicing. 

In addition to an implicit scheme, we also implemented an iterative explicit scheme, in which case the finite difference 
equations are written as 


a”+^ - a. 


^=<e”+/3r 


e”+^ - er 


X («+i/2)"«+i - <) - «-i/2)^« - <-i)) 


-ar Tarver 


Here the radius r denotes an areal radius in this Section (denoted with in the rest of the paper). Prior to excision, the 
boundary conditions (7-36) and (7-37) can be finite-differenced as 

-2-=-^ 
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and 


n+l 
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n+l 

2 


2 


V 


n+l 
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(1 + 2Mlr^/^Yn 


^n+1 
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n+l 
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Once the stellar radius is smaller than the excision radius i?ex we switch to the excision boundary conditions (7-38), which 
can be implemented as 


.,n+l 


, 11 - 1-1 


„n+l _ „n+l 


(B13) 

(B14) 


The outer boundaries are identical to those for Schwarzschild coordinates and maximal slicing. 

An iterative Crank-Nicholson scheme can be constructed in the following way (cf. Teukolsky 2000). In a predictor step, 
the above equations are solved for approximate new values 0 +^ and 6+^. The average between these values and the old 
values a" and e" is then used on the right hand side to compute improved new values 0 +^ and 6+^ in a first corrector 
step. These are again averaged with a" and e” and used on the right hand sides for the second and final corrector step, 
which yields the new values 0 +^ and 6+^. 

The Kerr-Schild results presented here were obtained with the above explicit scheme. However, we also implemented 
an implicit scheme very similar to that in Section B.l and obtained very similar results. In particular, we found the 
interesting result that excision boundary conditions (7-38) can be used with an implicit scheme without encountering 
problems. Implemented as in Section B.l, the implicit scheme is first order in time (but could easily be made second 
order), while the above iterative scheme is second order in time, which resulted in more accurate results in late-time 
evolutions. 
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